function [imageClass probClass] = classifyDigit(inputImage, atlas, normalized)
% classifyDigit: 
% Parameters: inputImage - the input image that we wish to classify
%             atlas - The cell containing the atlas data
%                altas.templateAll =  10*16*16 matrix which is the set of
%                10 image templates one for each class
%                atlas.cpAll = cell of size {10,1} with each member being a
%                matrix of size Ns*2*ncp
%                atlas.alphaAll = cell of size {10,1} with each member being a
%                matrix of size Ns*2*ncp
%             normalized - if "1" then account for class imbalance which we
%             have not thought about as yet :(
% Returns:  imageClass - The class the image belongs to

% for displaying debug and other info
verbose = 1;

% First register the input image with each template and find the momenta
% vectors. Store the deformation vectors in a 10*2*ncp matrix defMomAll.
% Index 1 contains the deformation vectors for class 1 while index 10
% contains the vectors for class 0.
defMomAll = zeros(size(atlas.alphaAll,1), size(atlas.alphaAll{1},2), size(atlas.alphaAll{1},3));
for i = 1:10
    source = squeeze(atlas.templateAll(i,:,:));
    % rescaling the image to lie between 0 and 1
    source = (source-min(min(source)))./(max(max(source))-min(min(source)));
    target = inputImage;
    [deformIgnore cpIgnore mom] = pairRegistration(source, target, squeeze(atlas.cpAll{i}(1,:,:)));
    defMomAll(i,:,:) = mom;
    if(verbose == 1)
        disp(['registered with class number: ',num2str(i)]);
    end
end

% With the momenta vectors compute the Mahalanobis distance of each of the
% deformations from the corresponding mean of the deformations of each
% class.
% Store it in a vector mDistFromClass of size 10*1 with mDistFromClass(1)
% being the distance of the image from mean of class of digit 1
mDistFromClass = zeros(10,1);
for i = 1:10
    % reshape the 2*ncp deformation vectors to a 1*(2ncp) vector
    defVect = [squeeze(defMomAll(i,1,:))' squeeze(defMomAll(i,2,:))'];
    % now reshape the reference data from Ns*2*ncp to Ns*(2ncp)
    clear atlasVect;
    atlasVect = zeros(size(atlas.alphaAll{i},1),size(atlas.alphaAll{i},2)*size(atlas.alphaAll{i},3));
    for j = 1:size(atlas.alphaAll{i},1)
        atlasVect(j,:) = [squeeze(atlas.alphaAll{i}(j,1,:))' squeeze(atlas.alphaAll{i}(j,2,:))'];
    end
    mDistFromClass(i,1) = mahalanobisDistNew(defVect, atlasVect, size(atlasVect,1));
    % normalize to account for imbalanced data
    if(normalized == 1)
        % not implemented / thought about!!
    end
end

mDistFromClass = 1./mDistFromClass;

% need to find the probability for each class
probClass = zeros(10,1);
for i=1:10    
	probClass(i,1) = mDistFromClass(i,1)/sum(mDistFromClass,1);
end

%find the image class which the inputImage is closest to
imageClass = find(probClass==max(probClass),1);

if(imageClass == 10)
    imageClass = 0;
end

end

function result = mahalanobisDistNew(observationData, referenceData, numRefSubjects)
% mahalanobisDistNew: This function finds the mahalanobis distance between
% the observation data which is a 1 x ndims matix and the reference data
% which is a matrix of dimensions nobs x ndims

% find the mean of the reference dataset
% subtract the mean from the datapoint in the observation dataset
uRef = mean(referenceData(1:numRefSubjects,:));
coVar = cov(referenceData(1:numRefSubjects,:));

result = sqrt((observationData - uRef)*pinv(coVar)*(observationData - uRef)');
   
end
